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We describe a new, microscopic model for diffusion that captures diffusion induced 
fluctuations at scales where the concept of concentration gives way to discrete par- 
ticles. We show that in the limit as the number of particles N — > oo, our model 
is equivalent to the classical stochastic diffusion equation (SDE). We test our new 
model and the SDE against Langevin dynamics in numerical simulations, and show 
that our model successfully reproduces the correct ensemble statistics, while the 
classical model fails. 
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I. INTRODUCTION 

Fluctuations in concentration become important when modeling systems with small par- 
ticle density, due either to small concentration or small spatial scale. They can have a signifi- 
cant effect on the average behavior of many diffusion-reaction systems. Fluctuations in chem- 
ical reaction kinetics result from both the intrinsic stochastic nature of chemical reactions 
and concentration fluctuations. Concentrations fluctuations, in turn, result from both reac- 
tion fluctuations and the the random, thermal motion of reaction species. Thus, reaction- 
diffusion systems couple kinetics and diffusion at both the deterministic (macroscopic) and 
stochastic (microscopic) scales. When nonlinear reactions are present, initial fluctuations 
can induce instabilities that lead to interesting macroscopic behavior such as pattern for- 
mation and oscillations [6J. When diffusion induced fluctuations are also present, we may 
see different behavior than deterministic models of the same systems predict [U [10], EH EH] • 
This highlights the importance of having good models for diffusion induced fluctuations. 

Various Lagrangian- and Eulerian-frame representations exist for theoretical and numer- 
ical modeling of reaction-diffusion systems at the fluctuation scale. One fundamental way 
to model reaction and diffusion is through Langevin dynamics, i.e. particle tracking. This 
is a Lagrangian-frame representation that tracks the motion of individual particles whose 
dynamics is described by the Langevin equation (possibly in over-damped form), and models 
reactions based on some probabilistic or deterministic function of inter-particle distance. 

In many situations it is more convenient to work in an Eulerian frame where one is in- 
terested in the concentration of material at a point in space (or the number of particle in a 
small volume). The classical Eulerian description of diffusion is the diffusion PDE which one 
can derive by considering an ensemble of particles in Brownian motion. However, this is a 
macroscopic model for average particle density, and does not include fluctuations. A meso- 
scopic description that includes fluctuations is the Multivariate Master Equation (MME), a 
spatially discrete (Eularian) continuous time Markov chain. The MME has been used to ob- 
tain some important rigorous results concerning the onset of instabilities in reaction-diffusion 
systems [13J. The advent of powerful computers has enabled numerical sampling of master 
equations to become feasible. Exact sampling methods, such as stochastic simulation algo- 
rithms (SSAs) exist, but are usually slow and, more importantly, progress in random time 
steps [9J. This complicates multi-scale modeling, especially where one would like to couple 
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microscopic to mesoscopic to macroscopic models where the transitions between mesoscopic 
and macroscopic regimes may change dynamically with space and time. 

One can derive a stochastic diffusion equation (SDE) from the MME as a thermodynamic 
limit. The SDE adds a stochastic flux to the classical diffusion equation. The SDE can be 
discretized and use in numerical simulations where one wished to model diffusion induced 
fluctuations. It also evolves in fixed time steps. Therefore, the SDE can seamlessly inte- 
grate with a forward- Euler finite-difference integration of the deterministic diffusion PDE. 
The SDE is little more expensive than for deterministic diffusion - generating Gaussian 
random variables being the additional expense. However, since the SDE represents the ther- 
modynamic limit, a valid theoretical question is to what degree models very small particle 
densities. 

We have found a new representation, called the multinomial diffusion equation, (MDE) 
that describes the evolution of the numbers of particles in a spatially discretized field in fixed 
time steps. Using numerical simulations, we compare the diffusion induced fluctuations in 
both the MDE and SDE to those observed in a particle tracking model. We find that our 
new MDE more closely reproduces diffusion induced fluctuations than the SDE. Therefore, 
we conclude that the MDE provides a theoretical middle ground between the MME and the 
SDE. We also found that the MDE can be used as an efficient and accurate finite difference 
method for modeling diffusion at the particle scale. It is comparable to a particle simulation 
in accuracy, yet is almost as computationally efficient as the SDE (in a finite-difference 
discretization), and also evolves (synchronously) in fixed time steps. 

For the remainder of this article, we will continue to use the term "diffusion" rather 
than "Brownian motion", which might more accurately specify that we are looking at the 
diffusion of particles, as opposed to heat, for instance. However, we will from now on use the 
term "particle density" instead of "concentration" to emphasize that we are in the regime 
of individual particles. 

II. MULTIVARIATE MASTER EQUATION 

The Multivariate Master Equation (MME) models the numbers of particles in M voxels 
of size Ax in a spatially discretized domain of size L = MAx [8]. The state of the system 
is a spatial field of particle numbers recorded in the vector N = [Ni,N 2 , ...,N M ] where iVj 



4 



is the number of particles in the i th voxel (centered at Ax{i + |)). In a transition event, a 
single particle hops from the j th voxel to the i th voxel. Such a transition changes the state 
from N to [JVi, N 2 , Nj - 1, N { + 1, iV M ]. More compactly, iV ->■ N + A, where A 
has only two nonzero elements: Aj = — 1 and A, = 1. Let V[A\N(t), At] be the probability 
that that the transition N(t) — >■ iV(t) + A occurs during the next small time increment At. 
We would like to have an expression for P(N(t + At),N(t)) = P[N(t + At)\N(t)] P(N(t)). 
Since N(t + At) = N(t) + A, P(N(t + At)\N(t)) = P(A\N(t + At)). Therefore, we have 



P(N(t + At),t + At\N(t)) (1) 

A 

i j 



where % — > j stands for Aj = — l,Aj = 1. As in pure Brownian motion, we will assume 
that individual particles do not interact. With this assumption, V[i — > j] depends only 
on Ni(t). Also, when the linear dimension of a voxel is larger than the mean free path 
of a Brownian particle, we need only include nearest neighbor hops, i.e. \i — j\ G 0,1 . 
The mean free path (A) for Brownian motion is a measure of the distance a particle can 
travel after an impulse from the surrounding fluid. A good estimate assumes the particle 
starts at thermal speed \ — — } giving A ~ ^gj where 1] is the fluid viscosity, and a is 
the particle diameter. A "particle" must be larger than molecular size ~ 10~ 10 m. Using 
molecular size, and the density of stone (~ 10 6 &;/m~ 3 ), the mean free path is ~ 10 _10 m 
- no more than the diameter of the particle itself! With these assumptions, and using 
V(i i) = 1 - V[i -)■ % - 1] - V[i -)■ % + 1] we transform Eq. Q to 
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p(,_l_ > ,| iV ,_ 1 (t))p(iV,_ 1 (t)) 
+ V(i + l^i\N i+1 (t))P(N i+1 (t)) 
+ (1-V(i^i- l\N t (t)) -V(i^i- l\Ni(t)P(Ni(t)) 



which has the form of a Chapman-Kolmogorov equation for the interval from t to t + At. 

An informal way to define a transition rate W[j — > i] from a transition probability 
V[j — > i) is to say that W[j —>■«] = jl'Pfj — >■ i]. One can make this rigorous when 
V[j — >■ i\ « gj^At + o(At) and P[i i] « 1 - <7i,iAt + o(At). In this manner, Eq. ([2]) gives 



P(N(t),t) 
dt 



£ (3) 

i 

W[j + l^j|iV i+1 (t)]P(iV i+1 (i) ) i) 
_ W [ i ^ i + l| iV . (t )]P( i V. (t ) ;t ) 

-W[j^j-l|iV,(t)]P(iV,(t),t) 
+ W[j - 1 J |iV J _ 1 (t)]P(iV^ 1 (t), t) 



which is known as the multivariate master equation for diffusion. 

The theory of continuous time Markov chains allows us to decompose this process into 
two independent random processes: (1) a random waiting time until the next transition and 
(2) a random selection of which transition occurs j5]. The distributions for these random 
numbers depend on the transition probabilities W(A\N). This is the basis of exact sampling 
algorithms such as the Gillespie algorithm [9]. 

III. MULTINOMIAL DIFFUSION EQUATION 

We now describe a representation that is spatially discrete (as is the MME), but evolves 
in fixed time steps (as does particle tracking). We use the same definitions as we did in 
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deriving the MME, except that we do not restrict to single particle exchanges. Instead, 
A = (Ax, A 2 , . . . , A M ) where Aj can have any value < Aj < N, so long as £V A, = 0. We 
again invoke the assumption that when Ax is larger than the mean free path, only consider 
nearest neighbor exchanges. Since we are now working with fixed time steps, we will use 
time steps as an index placed as a superscript. Let the vectors L* and R 1 record the random 
number of particles that jump out of voxel % to the left and right respectively (with suitable 
boundary conditions) in the time interval from t to t + At. 

Let nAt be the probability that an individual particle can jump into the next voxel 
during an interval of size At. In this case, the probability that L\ particles jump to the 
left out of voxel i, and R\ to the right is given by the multinomial multinomial distribution 
M(Nf,kAt,kAt): 

P\L* J*| = N ^At)^ikA^il-kAtfl-^ 

This is the essential feature of the MDE. 

Can we derive a master equation for the MDE? Let us define the shift operators L and 
R such that 



[LL], = [L} 1+1 (5a) 



= (5b) 



M]i = [RUi ( 5c ) 
\RR]i = [R]i-! (5d) 

— * 

and LM. = WL = I. For example, L pulls the value of L in slot i + 1 back to slot i, and 
likewise for R. This gives the relationship A = LL* + M.R t — (L* + R l \ Conditioning on L 

— * 

and R, and again using only nearest neighbor exchanges, we can write 
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P(N(t + At),t + At\N(t),t) = 

i 

+ V[Lj +1 ,Rl 1 \Nl 1 ]P(N t+1 ,t) 
+ {l-V[L\,R t i \Nt))P{Nt,t) 

(6) 

This also has the form of a Chapman-Kolmogorov equation, and is comparable to Eq. ([2]). 
However, from Eq. (|4]) we see that the transition probabilities in Eq. ^ are not o(At). 
Therefore, we can not construct a master equation as we did for the MME. 

Nevertheless, we still can generate exact realizations of N(t). At each time step, we can 
generate the random vectors Uit) and according to Eq. Q, and then perform the 

updates 

fit+At = fit + h £t + R pt _ + fit^ ^ 

We will also write this term-by-term 



Nl +At = (8) 
N- + L\ +l - R i>t - L\ t + R\_ x 

IV. STOCHASTIC DIFFUSION PDE 

There is a classical stochastic diffusion PDE that can be derived by various methods in 
the thermodynamic limit of iV — > oo. 

d -^dT = D &^ t] + Y x ^D^J)i^t) (9) 

Keizer derives Eq. (|9j) using thermodynamic potentials [12] . Gardiner derives Eq. ^ from 
MME using a Van Kampen system size expansion [8]. Ironically, these derivations employ 
the limit N — > oo, even though fluctuations are only significant when N ^ oo. This suggests 
there is a lower limit of particle density where this description will apply. For instance, we 
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wonder if Eq. (J9J) can accurately model diffusion induced fluctuations as well as a particle 
tracking model when the particle density is very small. 

V. THERMODYNAMIC LIMIT 

It is common to use a multidimensional Gaussian distribution to approximate a multi- 
nomial distribution p. A multinomial distribution P[ni,n 2 , . . . ,Um] = N\ Yli=i ttp w ith 
^2i=i n i = ^5 can be approximated by a multivariate Gaussian with means /ij = Npi, vari- 
ances S^j = Npi(l—pi), and covariances = —NpiPj. The approximation becomes better 
as A gets larger, but worse as each pj gets smaller. 

Let us consider what happens if we make this approximation in Eq. (J8l). Since k is a 
probability rate, as t — > 0, &Ai becomes very small. In this limit, £^ ps NkAt + o(At), and 
Sjj- o(At 2 ). The multinomial random variables in Eq. (j8j) become independent Gaussian 
random variables, and we have 



A* +At = A' + KAtlNtu - 2A* + JV*_J (10) 



From the deterministic part of Eq. (10) we learn that kAt = This is expected, since a 
particle has a high probability of traveling a distance Ax in a time interval At = ^ D j Ax 2 . 

Conservation of mass requires that we cannot remove more than Aj particles from voxel 
i in any time step. In the MDE, using the multinomial distribution ensures that mass 



is conserved. We ask how large must A be so that Eq. (10) will almost never violate 
conservation of mass? Conservation of mass for the deterministic part of Eq. ( 10 ) requires 

m 



(ii) 



9 



However, due to the fluctuating part of Eq. (1 10h there is a finite probability that jV* +A * < 



even if Eq. (11) holds. To obtain a rough estimate for how rare such an even would be, we 



require that the total number of particles that leave Ax/2 < x < Ax/2 in a very small time 
At is between and N by some number standard deviations, s. 



H - 2sa > 



N > 



As 2 

DAt 
Ax 2 



(12) 



We might expect that iV would increase quickly with s, as we require smaller and smaller 
probability of violating conservation of mass. However, in the continuum limit, At — > and 



Ax 0, thus Eq. (12) requires that — > regardless of how we take this limit. Hence, 



Eq. (12) shows that we have the more stringent requirement that N — > oo as At — > 0. Thus 



a continuous version of Eq. (|8j) is not valid for large N, but strictly for N — ?► oo. 

Leaving this matter aside for the moment, we will show how one can draw an equivalence 
between Eq. ^ and Eq. ^ using Eq. (10). To obtain an expression for particle density p 



rather than particle number N, we divide Eq. (10) by Ax, and obtain 



J+At 



Pi 

DAt, , 

Ax 2 Lr ' m 



[Pl+i ~ 2/4 



pU\ 



+ 



DAt 

Ax 3 



Pi+lCi+l V Pii\,l V Pi £i,2 + V Pi-l 



(13) 



Finally, we apply the identity cti £ + £ = a/°"i + °\ £ ^° combine some of the Gaussian 
random variables and we are left with: 



(14) 
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Taking the continuum limit of a discrete stochastic equation such as (14) is not trivial. 
Garcia et al. have derived a rigorous discretization of the SDE, Eq. Q [7\. We now consider 
an informal derivation of this same discretization. We start with the usual discretization for 
the deterministic part 



P ? At =pi+^ipii-2pi+pU] (15) 



The Ito time discretization requires a factor of v At for the fluctuating part 



P? At = Pl + ^lpli-2pl + Pli] (16) 



In the appendix, |XJ we show an informal way to discretize -§^y/f>V- Combining these 



parts, we obtain exactly Eq. (14) 



VI. SIMULATIONS 

The SDE, Eq. ([9]), is derived in the thermodynamic limit, i.e. where N is strictly infinite. 
On the other hand, the MDE is a particle scale model. To compare how well these two models 
reproduce diffusion induced fluctuations, we performed numerical simulations comparing 
the SDE and MDE to a particle tracking model, which we consider more fundamental and 
realistic. We generated realizations of diffusion for a total time T max in a periodic domain of 
length L, initialized with N particles distributed uniformly over the domain. For the SDE 
and MDE, we discretized the domain into N v voxels of size Ax = L/N v , creating an initial 
particle density of no = Nq/L. We also overlaid this grid on the particle simulation domain 
in order to calculate particle density. To study the effects of time step and particle density, 
we varied At and Nq while fixing D and Ax - which effectively defined our space and time 
units. 
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A. Particle Tracking 

In the particle tracking simulations we initially filled a domain of size L with N particles 
distributed uniformly. At each time step we used over-damped Langevin dynamics to update 
the positions of the particles 



< +At = x t n + V2DM£ t n +At n = l...N (17) 

At each time step we counted the number of particles in each of the N v voxels defined above 
to obtain the particle density. 



B. MDE 



We evolve the MDE using Eq. (J8J). Realizations of the MDE require generating multino- 
mial random variables. To generate multinomial random variables ni and n 2 from N with 
probabilities q\ and q 2 , we used a sequential approach based on successive binomial random 
variables. To generate two multinomial random variables from the multinomial distribution 
Ai[N, qi, q 2 ], we first chose ri\ from from the binomial distribution B[N, qi], and then chose 
n 2 from the binomial distribution B[N — nx, 92/(1 — qi)}. 



SDE 



The discretization in Eq. ( 14 ) has a finite probability of generating negative concentra- 
tions. Taking the absolute value of the concentration would clearly create a bias in the mean. 
To minimize the bias, we allowed the concentrations to be negative in the deterministic part, 
but took the absolute value for the square root in the fluctuations. However, this approach 
would not work in simulations with reactions. Allowing negative concentration would prop- 
agate through any reaction channel with an odd order, possibly leading to runaway negative 
concentrations. Using this correction, we used the following discretization to integrate the 
SDE 
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Pi 



t+At 



f DAt r . , , n 

Pi + -^[pt+i-M + Pi-i] 



(18) 



DAt 

~Kx 2 



Pi+l+Pl\Q,l- \^\Pi + Pi-x\Q,2 



where x and £* l are two different IID Gaussian random variables generated for voxel i at 
time t. 



VII. RESULTS 

The statistical properties of our models are seen in the ensemble statistics. We will use 
the notation for ensemble averages (over realizations u) and for spatial averages 
(over L). Let n a (xi) be the particle density at point Xi in the a th realization uj a . We define 
the ensemble mean particle density as fi(xi) = (na^Xi))^, and ensemble particle density 
fluctuations as cr 2 (xi) = ((n a (xi) — /u(xj)) 2 ) w . We also define a more concise number, which 
is the spatial domain averaged fluctuation a 2 = (a 2 (xi)) n . fi(xi) should match the analytical 
steady-state solution of the deterministic diffusion equation. However, we have no analytic 
expression for a 2 (xi). Of our three simulation methods, the particle method is the most 
fundamental and realistic. So we measure the "convergence" of fi(xi) and a 2 in the MDE 
and SDE simulations by how well their statistics match those generated in the particle 
simulations. 

Fig. [T] shows a typical result in which we see the fluctuations in the MDE simulation 
being somewhat larger than in the particle simulation, and the fluctuations in the SDE being 
larger still. The legend box gives ((/x(x) — no) 2 ) n as a measure of how well the the mean 
ensemble means agree - a test of convergence to the analytical steady-state solution. Fig. [2] 
shows a 2 for different values of At and uq. This plot clearly shows that Aa 2 /A(At) — > as 
At — > 0, however only the MDE converges to the same value of a 2 as the particle method. In 
other words, the MDE closely replicates the ensemble statistics of a particle simulation. On 



the other hand, to the extent that a numerical integration of Eq. ( 14 ) represents a solution 
ofEq. @, our data also suggests that the SDE, Eq. ([9]), equation is not an adequate model 
for very small particle density. In fig. |3j we see that the accuracy of the SDE appears to 
break down at about one particle/voxel. We suspect this breakdown will occur at higher 
particle density in higher dimensions. 
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FIG. 1. Typical results for an ensemble (N = 8192) rangeof simulations for DAt = 0.25, no = 0.5, 



Tmax = 32. The legend gives a measure of fit to analytical solution (see text). 



VIII. DISCUSSION 



From our simulation algorithms, we see that generating realizations of the MDE operates 
in a very similar way to typical finite difference methods used to solve the diffusion PDE 
or SDE. This suggests possible applications of this model for numerical simulations where 
one needs to model some spatial regions at the particle scale and include diffusion induced 
fluctuations. 

The grid-based MME is a somewhat more fundamental model than the MDE, and being 
expressed as a master equation may make it more amenable to some analytical work. From 
a practical point of view, however, the MME has some limitations. Exact sampling meth- 
ods, such as the Gillespie algorithm, evolve in single particle events with with extremely 
small, random time steps. Suppose a multiscale simulation has two or more disjoint regions 
requiring particle scale resolution. Using the MME, each of these regions will produce in- 
dependent, tiny time steps. However, finite difference, smoothed particle hydrodynamics, 
finite element, and most other prevalent techniques for modeling spatiotemporal fields op- 
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FIG. 2. Convergence of ensemble fluctuations as a function of time step (At). For large particle 
density, all models converge to the same value as the particle model as time steps decrease. However 
for small particle density, the Gaussian distribution models do converge, but to the wrong value. 

erate in synchronous time steps. Not only will the MME regions become the computational 
bottleneck, it will be difficult to couple the independently asynchronous MME regions to 
synchronous regions. Also, we previously discussed how the issue of negative concentrations 
causes trouble for the SDE. The MDE, on the other hand, can integrate well with other 
grid-based models that operate in synchronous time steps, and requires no corrections for 
negative concentrations. 

Recently, Alexander et al. proposed using Eq. ^ to couple between particle-based grid- 
based methods for simulating diffusion [3]. However, if one could approximate a particle- 
based simulation with a much less computationally expensive grid-based method, this could 
be even more valuable. With this suggestion, we should at least minimally address the 
issue of performance. For a simulation with M particles, the computational expense of one 
time step is roughly M. For a simulation with N grid points, the computational expense is 
roughly aN, where a ~ 2 d is some geometrical factor depending on the spatial dimension d 
and the type of discretization. Considering the particle density n = M/N, we see that a 
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FIG. 3. At the smallest value of At we studied, all methods converge to the same ensemble 
fluctuations for large particle density, but diverge as smaller values. 

grid-based simulation is less expensive when M > aN, i.e. Uq = M/N > a > 2. In this case, 
the critical particle density where the MDE has a performance advantage is independent of 
both M and N. 

In a complex system one may not know, before hand, what the particle density will be 
at every point in the domain at every time. Using the MDE spares one from having to 
dynamically create and link particle regions to grid regions, constantly performing checks 
to see where and when and where to put the interfaces. More importantly, if one includes 
reactions, the performance advantage increases. With a bimolecular reaction, a particle sim- 
ulation requires an additional ~ M 2 operations per time step; a grid-based simulation incurs 
only another ~ N operations. Here, a grid based simulation such as the MDE outperforms 
a particle simulation as long as M 2 + M > N(a + 1). When many particles are involved 
in the simulation, M 2 ^> M, and (roug hly) M > y/(a + l)N, giving n > y/{a + 1)/N. In 
this case, the critical particle density above which the MDE gains a performance advantage 
shrinks very steeply with the number of voxels. For even a modest number of voxels, this 
quickly approaches n < 1, the regime where the MDE becomes superior to the SDE in 
accuracy. 
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Monte-Carlo simulations such as the ones we have discussed are intended to generate 
actual realizations of a physical model. To insure this, we must adhere to proper space and 



time ordering. This means that E(x,x + Ax) in Eq. ([8]), and \/p + + p£, in Eq. (14) must 
be the same number when used for updating the particle density at points x and x + Ax. 



Comparing Eqs. ( 14 ) and (|9j) we see that proper space ordering also relates to conservation of 
mass: the number of particles crossing between x and x+Ax is the same when looked at from 
either side. Furthermore updating neighboring voxels with two different and independent 
random numbers would lead to enhanced flux fluctuations. One often wishes to accelerate a 
finite difference simulation with implicit or semi-implicit methods. It might be worthwhile 
to investigate if proper space and time ordering limits this. 



IX. CONCLUSION 



We have proposed a new discrete model of diffusion (the MDE) that sits between the 
multivariate master equation (MME) (a mesoscopic model) and the the classical diffusion 
SDE (a field model) - a microscopic, particle model occupying an even smaller regime. Using 
this model we give another, perhaps more intuitive, derivation of the classical diffusion SDE. 
We perform simulations showing that at very small particle densities, the MDE very closely 
approximates the statistical properties of a particle simulation, while the SDE does not. To 
our knowledge, this work represents the first time the classical diffusion SDE has been put 
to the test against a particle model. 

In addition to a new theoretical model for diffusion, we suggest important practical 
applications of the MDE. Although the MME is a lower-level description, it is not well 
suited for being coupled to other grid-based methods. The MME is often very slow, and each 
domain being modeled by the MME would generate its own time step. Particle simulations 
can be more efficient for modeling diffusion for small particle densities when no reactions are 
involved, but do not integrate seamlessly with finite-difference methods. When even simple 
reactions are involved, particle methods become much less efficient, than the MDE. 



Appendix A: Discretization of ^-i//^ 

Our goal is to find a numerical discretization for 
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^y/p(x,t)ri(x,t) (Al) 

where rj(x) is unit, Gaussian, white noise. 

An accepted way to form a decent discretization for the derivative of a function / = dF / dx 
is to use a centered difference. One can think of the centered difference discretization as 
estimating the slope by taking the difference of average values of the function to the left 
and right of a point: 

,n F(x+Ax)+F(x) _ F(x)+F(x-Ax) 

4~ ~ 2 (A2) 

dx Ax v ' 



We would like to use this approach to discretize Eq. (Al). However, due to the noisy 



fluctuations, Eq. (A2) may not be a good enough estimate of the average. To be more 



confident, we could use 



dF £ C ^ F Wx' - £ j:_ Ax F(x>)dx> 



dx 



2Ax 



(A3) 



i px+Ax i i-x 

2Ax / F{X ' )dX ' ~ 2Ax / a F{X ' )dX ' 

J x Z,L-\X J x -Ax 



In our case, F = y/prj. To use Eq. (A3), we need to know how to integrate yfpq. For 
this we turn to a method due to Chandrasekhar. Chandrasekhar approximates the integral 

j h a f(x)rj(x)dx by 



f(x)rj(x)dx 



N 



\ i=i 



(A4) 



where Ax = and Xi = a + Ax(i — 1). We now approximate the first integral in Eq. (A3 ) 



by setting N = 2, a = x, and b = x + Ax, and the second integral by setting a = x — Ax, 
and b = x. In principle, it makes sense to index r\ at any point between x and x + Ax. In 
this way, we finally obtain 
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